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ABSTRACT 

An explanation of long-lived Saturn’s North Pole hexagonal circumpolar jet 
in terms of instability of the coupled system polar vortex - circumpolar jet 
is proposed in the framework of the rotating shallow water model, where 
scarcely known vertical structure of the Saturn’s atmosphere is averaged out. 
The absence of a hexagonal structure at the Saturn’s South Pole is explained 
along the same lines. By using the latest state-of-the-art observed winds in 
Saturn’s polar regions a detailed linear stability analysis of the circumpolar 
jet is performed (i) excluding (“jet-only” configuration), and (2) including 
(“jet+vortex” configuration) the north polar vortex in the system. A domain 
of parameters: latitude of the circumpolar jet and curvature of its azimuthal 
velocity profile, where the most unstable mode of the system has azimuthal 
wavenumber 6, is identified. Fully nonlinear simulations are then performed, 
initialized either with the most unstable mode of small amplitude, or with the 
random combination of unstable modes. It is shown that developing barotropic 
instability of the “jet+vortex” system produces a long-living structure akin to 
the observed hexagon, which is not the case of the “jet-only” system, which 
was studied in this context in a number of papers in literature. The north 
polar vortex, thus, plays a decisive dynamical role. The influence of moist 
convection, which was recently suggested to be at the origin of Saturn’s north 
polar vortex system in the literature, is investigated in the framework of the 
model and does not alter the conclusions. 
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1 Introduction 


Visible imagery acquired during the Voyager 1 & 2 fly-by of Saturn in 1980-1981 
revealed the hexagonal cloud pattern of Saturn’s north pole at about 77°N plane- 
tographic latitude (Godfrey, 1988). As followed from these early observations, the 
motion of the associated cloud structures suggested that the structure is related to 
a circumpolar jet-stream, but the hexagon pattern itself appeared to be stationary 
(i.e. very close to the rotation of Saturn’s interior). The early diagnostics were con- 
firmed first in the 1990s by Hubble Space Telescope observations (Sanchez-Lavega 
et ah, 1993), and then in the 2000s by the Cassini orbiter firstly in thermal infrared 
images (Baines et ah, 2009), followed by visible imagery once Saturn’s northern 
pole came out of the winter darkness. The high spatial resolution of the Cassini 
visible and infrared images permitted to evidence, in addition to the hexagonal 
jet, the presence of the north polar vortex (NPV) . The exceptional duration of 
the Cassini mission allowed Sanchez-Lavega et al. (2014) to conclude that the 
hexagon resists seasonal changes, before Antunano et al. (2015) compiled the re- 
peated cloud-layer observations of the hexagon to obtain a complete profiling of 
the winds in this structure. Thus far, such hexagonal feature has not been observed 
at Saturn’s South pole. 

The persisting hexagonal pattern in polar projected maps demonstrates that a 
prominent wavenumber 6 perturbation shapes Saturn’s circumpolar jet stream at 
latitude ~77°N. An early interpretation by Allison et al. (1990) was that the anti- 
cyclonic North Polar Spot (NPS) vortex, visible at that time in the vicinity of the 
hexagon forced a stationary planetary (Rossby) wave which gave the polar jet its 
hexagonal shape, but the apparent disappearance of the NPS between 1995 and 
2004 invalidated this scenario (see Sayanagi et al. (2016) for the history of the 
observations of the NPS and the hexagon). The nature and origin of the hexagon 
has been addressed since then both in laboratory experiments and numerical mod- 
els. Polygonal structures are often reported in rotating flow experiments, although 
the wavenumber 6 configuration remained elusive (Vatistas et al., 1994; Marcus 
and Lee, 1998; Jansson et al., 2006; Bergmann et al., 2011). Aguiar et al. (2010) 
discussed more specifically the relevance of rotating tank experiments to modeling 
Saturn’s polar jet and concluded that the hexagonal structure could be resulting 
from the barotropic instability of the jet, and that the wavenumber 6 prominence 
is conditioned by the intensity of the jet and bottom friction. They were first to 
notice a dependence of the azimuthal wavenumber of the most unstable mode on 
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the deformation radius in their instability analysis in the framework of the quasi- 
geostrophic (QG) model. Using a general circulation model with a domain limited 
to Saturn’s northern polar region, and a polar jet as initial condition, Morales- 
Juberias et al. (2011) reproduced the structures found by Aguiar et al. (2010) in 
laboratory experiments: the hexagonal shape resulting from a vortex street formed 
by developing barotropic instability of the jet, with cyclonic and anticyclonic vor- 
tices at poleward and equator- ward sides of the jet, respectively. Yet, the vortices 
forming the street were too large and too strong, with respect to observations, and 
this mechanism was giving rise to larger propagation speed than the observed one 
(Sanchez-Lavega et ah, 2014). (Morales- Juberias et al., 2015) proposed an alter- 
native “meandering jet” model which matches the morphology (including sharp 
potential vorticity, hereafter PV, fronts) and phase speed of Saturn’s hexagon, 
provided that an ad hoc vertical shear of the jet is introduced. The same authors 
concluded that deep jets evolve into vortex streets and shallow jets evolve into 
meanders, both being possibly at the origin of the absence of seasonal variability 
of the hexagonal jet. Although the meandering jet model of Morales-Juberias et al. 
(2015) offers thus far the best agreement with the characteristics of the observed 
Saturn’s hexagon, the putative meridional temperature gradients accompanying 
the vertical shear of the (supposedly shallow) polar jet are yet to be confirmed. 
Despite these developments, there is not enough convincing dynamical explanation 
of the existence and origin of Saturn’s North pole hexagon, as well as the absence 
of its counterpart at Saturn’s South pole. In the present paper we propose an al- 
ternative approach to the problem. Our analysis is performed in the framework of 
a simple rotating shallow water (RSW) model, which is of use for modeling atmo- 
spheres of giant planets, cf. (Dowling, 1995). A specificity of our approach is that 
the model is considered in the polar tangent plane approximation, the so-called 
gamma-plane, in order to take into account the variations of the Coriolis parameter 
with latitude in the polar regions. The use of the model, which may be obtained 
by vertical averaging of the full primitive equations, is justified by the fact that 
information on the vertical structure of Saturn’s atmosphere is scarce. All vertical 
structure being averaged out, the model is barotropic, in the sense that it does not 
contain any vertical shear, although it does allow for vertical stretching of fluid 
columns. The model combines simplicity with dynamical consistency and, with 
the astronomical parameters of Saturn being given, contains a single adjustable 
parameter: the effective Rossby deformation radius, Rj- The QG model used for 
theoretical analysis in Aguiar et al. (2010) can be obtained from RSW in the limit 
of vanishing Rossby numbers. The RSW model, which allows for efficient high- 
resolution numerical implementations, has been used for the analysis of various 
features of the Saturn’s atmosphere: mid-latitude jets (Showman, 2007), equato- 
rial super-rotating jet (Scott and Polvani, 2008), NPV (O’Neill et al., 2016) (in 
two-layer version with a moist-convective forcing in the latter work). Yet, to our 
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knowledge, there are no studies in the literature applying the RSW model to the 
developing barotropic instability of the circumpolar Saturn’s jet. The particularity 
of our approach is that we explore the dynamics in two different configurations 
of the mean zonal velocity profile: (i) considering only the 77°N circumpolar jet, 
with observed velocity profile, and no central vortex at 90°N (i.e. the traditional 
“jet-only” configuration considered in all above-cited papers), and (ii) considering 
the full “jet + vortex” system with the observed zonal velocity profile, which was 
not considered in the existing literature in this context. The state-of-the-art mean 
zonal velocity profiles on Saturn (Antunano et ah, 2015) are used in linear stability 
studies and for initializations of nonlinear simulations. 

The paper is organized as follows: in section 2 we introduce the model, discuss its 
general features and vortex solutions, present the analytical fits of the observed 
velocity profiles, and formulate the linear stability problem. Section 3 is devoted to 
the study of the jet-only configuration, with presentation of the results of the linear 
stability analysis and nonlinear simulations of the saturation of the instability. In 
section 4 we analyse the jet+vortex configuration along the same lines. Finally, in 
section 5 we present our conclusions and discussion, including the influence of the 
effects of moist convection upon the evolution of the instability. 


2 Rotating shallow water model for large-scale atmospheric motions, 
vortex solutions and their (in)stability 


2.1 The model 


Rotating shallow water equations, in the absence of forcing and dissipation read: 

7^7 + (/ + ~)(z x v) = -gVh, (2.1) 

L>t r 

Dh 

=^ + W.v = 0. (2.2) 

Here, by anticipating the application to polar jets and vortices, we use the cylin- 
drical coordinates on the plane (r, 9), and z is the unit vector normal to the plane. 
v = ur + v6 is the horizontal velocity with u, v denoting radial and azimuthal 
components, respectively, D/Dt = d/Ot + v • V denotes the material derivative, 
V = fd r + 6{l/r)dg is the gradient in the plane, where d r and dg denote partial 
derivatives with respect to corresponding arguments, h is the thickness of the layer, 
g is gravity and / is the Coriolis parameter. 
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Bottom topography b(r, 9) can be easily introduced in the model by replacing h by 
h — b in (2.2), and can be used to parameterize the deep layers of the atmosphere. 
This would however introduce ad hoc parameters, which we want to avoid. 

The RSW equations possess a Lagrangian invariant, the potential vorticity (PV): 

5 • (V x v) + f 
Q= h ’ 

where z ■ (V x v) is relative vorticity. The PV anomaly 

PVA = q-f 0 /H 0 , (2.4) 

where fo and Hq are reference values of Coriolis parameter and thickness of the 
layer (see below), will be used for diagnostics in what follows. We consider the 
RSW equations in the polar tangent plane, the so-called q-plane, where the effects 
of sphericity in the Coriolis parameter are taken into account in the lowest-order 
approximation: 

f = fo sin(<#>) = /|| (l - + ■■■) = fo (l - 7‘f* 2 + •••) ; 0<r*<f. 

(2.5) 

Here / 0 = 2 O s , cj) is planetographic latitude, fl s and R s represent Saturn’s angular 
velocity and radius, respectively (we will neglect the effects of non-sphericity in 
what follows). Approximate values for the mean Saturn radius and fo are respec- 
tively 55000 km and 3.2 x 10 -4 s _1 (Read et al., 2009b; Sanchez-Lavega et ah, 
2014). As is well-known, the rotating shallow water system possesses an internal 
scale, the Rossby deformation radius Rj = sfgHo / f 0 . 

Equations (2.1), (2.2) can be obtained by vertical averaging between two material 
surfaces of the full three-dimensional primitive equations for the atmosphere in 
pseudo- height isobaric coordinates, e.g. (Bouchut et al., 2009), when one of these 
surfaces is considered fixed (a constant pressure level), and one is free. The poten- 
tial temperature is considered uniform through the layer in this approximation, 
although an extension of the model including horizontal potential temperature 
gradients is possible. In this case g is the gravity acceleration of Saturn, and H 0 is 
the thickness of the layer. The model can be also considered as a two-layer model 
with strong disparity in depths and/or densities of the layers. In this case H 0 plays 
a role of equivalent depth, and g is so-called reduced gravity, i.e. gravity weighted 
with the stratification parameter. The deformation radius in this interpretation 
is related to stratification properties of the atmosphere. Whatever the interpreta- 
tion is, at a given rotation rate /o/2 , the deformation radius Rd is the only free 
parameter of the model, as the velocity scale can be taken to be gH () . 


(2.3) 
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In terms of velocity components, the momentum and mass conservation equations 
of the model read: 


- — - Jv = -gd r h, 

Dt r 

(2.6) 

Dv uv 

(2.7) 

jrr + — +ju = -gd e h, 

Dt r 

d t h + - d r (hru ) + -dg(hv) = 0. 

(2.8) 


2.2 Vortex solutions 


Stationary axisymmetric vortex solutions of the equations (2.6) - (2.8) are those 
with zero radial velocity u = 0, some axisymmetric distribution of azimuthal ve- 
locity v = V(r), and the thickness field H{r ) related to V(r) through the cyclo- 
geostrophic (gradient wind) balance: 

— + fV = gd r H. (2.9) 

r 

For a given V(r), corresponding thickness profile H(r ) can be obtained by inte- 
gration of equation (2.9). Note that in the QG model used for theoretical analysis 
in Aguiar et al. (2010) the centrifugal acceleration (the first term in (2.9)) is ne- 
glected. While this approximation is fully justified for the circumpolar jet, it is 
less so for the central vortex. Herein we consider solutions of (2.9) with V(r) cor- 
responding to the observed time- and zonally averaged profiles of Saturn’s zonal 
wind (Antunano et ah, 2015). These profiles, with error margins, are represented 
by dashed lines in Figure 1. 

It is useful to have a simple analytic expression mimicking the observed velocity 
profiles. Although the linear stability analysis can be accomplished directly with 
digitized observational velocity profile, it is convenient to control the shape of 
the velocity profile in terms of a small number of paramaters. Dependence on 
these parameters could be then analysed. The PV, which is a key characteristic of 
the flow, is very noisy, if determined directly from the observed profiles, cf. Fig. 
2, whereas a PV field computed from the analytic fit enables a straightforward 
analysis of possible instability, see below. For initialization of the non-linear RSW 
simulations, the use of the analytic fit allows for straightforward integration of 
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Figure 1. Observed zonal velocity profiles in the northern (left panel) and southern 
(right panel) hemispheres and their fits used in the paper. Three dashed lines represent 
the observed non-dimensional mean and error margins of the averaged zonal velocity 
(Antunano et al., 2015), and the solid line is the composite of two model profiles (in 
red) fitting the NPV and the circumpolar jet. The fit of the velocity profile of the 
NPV gives e — 0.15, a — 0.42, fd — 1.3, ro = 0,m = 1 and that of the circumpolar jet 
gives e = 0.08, a — 0,/3 = 2,ro = 3.37, m — 3. Similarly, south polar vortex corre- 
spond to e = 0.16, a = 0.5 , (3 = 1.2, ro = 0,m = 1 and southern circumpolar jet jet to 
e — 0.08, a — f), fd — 2,r^ — 4.6, rn = 3. R ( i ^ 3200 km is taken for both poles. 




r r 

Figure 2. PV distributions corresponding to observed velocity profiles of Fig. 1 in the 
northern (left panel) and southern (right panel) hemispheres and their analytic fits. 

(2.9), and thus diminish the discretisation errors. Therefore, we fit the peaks of 
the velocity distribution with the help of the following simple formula: 

V (r) = e(r — ro) a e~ m ^ r ~ r °' >13 , a,/3,e,m>0, (2-10) 

where e measures the intensity of the velocity field, ro tunes the distance of the 
velocity peak of the jet from the pole, and other parameters allow to fit the shape 
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of the distribution (r o = 0 for the central vortex and a = 0 for the jet). All velocity 
profiles are normalized by their maximum values, so the maximum value of velocity 
is equal to e. By applying this formula to both central vortex and circumpolar jet 
(red curves in Fig. 1), we get synthetic profiles represented by solid black lines in 
the figure. We do not seek to reproduce the “shoulder” in the velocity profile of the 
central vortex, as it turns to be inessential for the dominant long-wave instability 
of the system which we are interested in (we made a test study to check this). 

The stability of any velocity profile V (r) (and corresponding H (r) obtained from 
the cyclo-geostrophic balance) with respect to small perturbations can be analyzed 
by standard means. We represent each field as the solution in question and a small 
perturbation (denoted by a prime): 

u(r,9,t) = u'(r,9,t), (2-11) 


v(r, 9,t) — V (r) + v\r, 9, t ), 


( 2 . 12 ) 


h(r,9,t) = H(r) + h'(r,9,t), (2-13) 

and inject these expressions in the full system (2.6) - (2.8). Retaining only the linear 
contributions in perturbation terms, we are looking for normal-mode solutions with 
harmonic dependence on time and polar angle 

(u' ,v' ,h')(r,9,t) = Re[(iu, v, h)(r)^ l9 ~ ut \ (2-14) 

where / and co are the azimuthal wavenumber and the complex eigenfrequency, 
co = ujr + icoj. We thus formulate the problem as an eigenvalue problem for eigen- 
frequencies co. A positive imaginary part for an eigenfrequency coj corresponds to 
an instability with a linear growth rate a = coj. 


Following common practice, it is convenient to formulate this eigenproblem in non- 
dimensional terms. By using the deformation radius Rd as the scale for r, f ( J 1 as 
the scale for t, and \JgHu as the scale for velocity, the eigenproblem in question 
takes the following form: 


IV* 

(f + ^r + Or-v) 

«'/),- + Vl,-(r- H-) 



IH* 




-D r * 

IV* 


1 

r- 
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u 


u 

X 

V 

= (jj 

V 


fl 


f] 


(2.15) 


where variables marked with the asterisks are non-dimensional, and D r * denotes 



the differentiation with respect to r*. The non-dimensional 7* = 0(10 3 ). It should 
be noted that the parameter e acquires in this way a meaning of the Rossby number. 

Pseudospectral collocation method (Trefethen, 2000) is used for linear stability 
analysis of the eigen- problem (2.15). The system is discretized over Appoint grid 
and D r * becomes the Chebyshev differentiation operator. To avoid the Runge 
phenomenon, Chebyshev collocation points are used, following Lahaye and Zeitlin 
(2015) where the method of Boyd (1987) was adapted to the stability problem of 
circular vortices. 


3 Instability of the “jet-only” configuration and its nonlinear satura- 
tion 


In this section, we analyse the instability of the North Pole circumpolar jet in 
the absence of the central vortex (“jet-only” configuration), following the existing 
studies in the literature (Aguiar et al., 2010; Morales- Juberias et al., 2011, 2015), 
albeit in a different model. The velocity profile of the jet is fitted by a Gaussian, 
i.e. the profile (2.10) with a = 0. The non-dimensional amplitude of the jet is 
e = 0.08. The values of parameters m and r 0 corresponding to the most reasonable 
fit of the observations with error bars (see Fig. 1) are m = 3 and r 0 = 3.37. 


3. 1 Results of the linear stability analysis 


We present in this subsection the results of linear stability analysis for the “jet- 
only” configuration, carried out along the lines of section 2. The barotropic in- 
stability of zonal jets, as well as its nonlinear saturation, are well understood in 
geophysical fluid dynamics. In the framework of the RSW model sufficient condi- 
tions of stability of zonal jets are given by Ripa’s criteria (Ripa, 1990). Like the 
famous Rayleigh - Kuo stability criterion for the quasigeostrophic flows (which 
should be adapted, however, to the current 7- plane context) the meaning of these 
criteria is that in order to have the instability, the PV gradients should change 
sign within the flow. Physically, this means that there should be a possibility for 
Rossby waves, which owe their existence to PV gradients, to propagate in opposite 
directions in the flow, phase-lock, and grow in amplitude. As follows from Fig. 2 
there is indeed a change of sign of the gradient of PV, which is by the way more 
pronounced in the Northern hemisphere, and hence we expect the barotropic in- 
stability. The main question then is: under what precise circumstances the mode 
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with azimuthal wavenumber 1 = 6 is the dominant instability mode of the circum- 
polar jet? Our linear stability analysis demonstrates that for the location of the jet 
at the distance from the pole r 0 corresponding to the observations, the dominant 
instability mode has an az im uthal wavenumber 1 = 6 for a choice of R ( i ~ 3200 km. 
This value of Rd is in line with the values of 2500—3000 km adopted in Aguiar et al. 
(2010) and O’Neill et al. (2016) (cf. their appendix 2), and are of the same order 
of magnitude as the estimates from observations in Read et al. (2009a), although 
about twice larger. Yet Rd remains to be better constrained by observations. 

Changing the radial (latitudinal) position of the jet, r 0 , changes the wavenumber of 
the most unstable mode. This is shown in Fig. 3, where the growth rates of unstable 
modes as a function of e and vq are displayed in a range of values of e and tq covered 
by our linear stability analysis). Fig. 3 shows that the azimuthal wavenumber of 
the most unstable mode varies from l = 2 to l = 8 with r 0 increasing from 2 to 5 
in non-dimensional terms. For some values of r 0 , several modes exhibit very close 
growth rates, especially at small e: in such configuration a combination of these 
modes would eventually shape a circumpolar jet without pronounced polygonal 
pattern. It is worth noting that an increase of velocity amplitude of the jet, e, does 
not change the azimuthal wavenumber of the most unstable mode, but increases 
the growth rate, cf. Fig. 3. 

It should be emphasized, however, that even a small difference in growth rates 
yields exponentially- growing differences between corresponding wave patterns, whence 
the primary importance of the leading mode. 

In Figure 4 we present the regions in the parameter plane tq — m where unstable 
modes with a given l have the highest growth rate at fixed Rd . Higher values of m 
correspond to stronger curvatures of the velocity distribution. The value of Rd also 
influences the results. For larger values of Rd, the observed position of the maxi- 
mum jet velocity shifts closer to the pole in non-dimensional terms and, as follows 
from the results described in the previous paragraph, the azimuthal wavenumber of 
the leading instability becomes smaller, and vice-versa. By reversing the argument, 
from our analysis and from the fact that the observed dominant wavenumber is 
l = 6, we can infer an estimate for Rd ~ 3200 km. 

An important unsolved problem is the difference in morphology between the hexag- 
onal northern polar jet and the southern polar jet which does not exhibit a pro- 
nounced polygonal form in Saturn’s atmosphere, as follows from the observations 
(Antunano et al., 2015). Our linear stability analysis of “jet-only” configuration 
provides a clue for explanation. As follows from Fig. 1, at a given value of Rd, the 
value of ro is higher in the southern hemisphere (i.e. the southern jet is farther 
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Figure 3. Dependence of the growth rate of the unstable modes on azimuthal 
wavenumber for different latitudinal positions of the jet. V(r) is Gaussian, with 
m = 3, a = 0, /3 — 2, « 3200 km. 

from the pole), which yields a leading instability with l > 8 (Fig. 3), and hence no 
hexagonal pattern at Saturn’s South Pole, although a partial hexagonal structure 
was discussed in the literature (Vasavada et ah, 2006). 


3.2 Non-linear saturation of the instability 


In order to study the nonlinear saturation of the instability identified in the pre- 
vious subsection, we perform direct numerical simulations with a well-balanced 
entropy-satisfying finite-volume scheme (Bouchut, 2007) for the RSW equations. 
The numerical scheme was successfully tested in studies of the saturation of the 
barotropic instability of jets and vortices (Lambaerts et ah, 2011; Lahaye and 
Zeitlin, 2015). The scheme has no explicit dissipation, only a numerical one. The 
dissipation is concentrated in the zones of strong velocity and thickness gradients. 
For flows with low Rossby numbers, as is the case here, the scheme is practically 
dissipationless and allows for high-resolution long-time simulations. The simula- 
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Figure 4. Left panel: Distribution of azimuthal wavenumbers of the most unstable modes 
in the plane of parameters r$ — m. Black (white) square corresponds to the observed pro- 
file in the northern ( southern ) hemisphere. V (r) is Gaussian, with e = 0.08, a = 0, (3 = 2. 
Azimuthal wave-numbers larger than eight are not included. Right panel: Dependence 
of the growth rates of the unstable modes of the “jet-only” configuration in northern 
hemisphere of the value of equivalent depth. Variation of the growth rate with azimuthal 
wavenumber for the jet corresponding to the black square in the left panel corresponds 
to Rd ~ 3200 km. 

tions below are initialized with the most unstable mode of small amplitude (several 
per cent with respect to the background flow) superimposed onto the jet, or with 
an ensemble of unstable modes with random phases. 


The evolution of the instability follows the well-known scenario of the saturation 
of the barotropic instability. Namely, a series of vortices of opposite sign appear 
at the crests and lows of the initial unstable wave, cf. Fig. 5. This is similar to the 
vortex-street patterns observed in simulations of Morales- Juberias et al. (2011) and 
laboratory experiments of Aguiar et al. (2010). As time goes on, the vortices inten- 
sify, resulting in a strong deformation of the polar jet. Although the outcome does 
preserve a six- fold symmetry of the initial unstable mode, the observed hexagon 
with quasi-straight sides (corresponding to sharp PV gradients) is not reproduced. 
Moreover, as follows from Fig. 6, the azimuthally averaged velocity of the result- 
ing structure displays formation of counter-jets at both sides of the original jet, 
which is typical for the saturation of the barotropic instability (e.g. (Lambaerts 
et al., 2011)), and does not match the observations. At later stages the hexagonal 
structure practically disappears (cf Fig. 7). The evolution of radial distribution 
of zonally averaged PVA, as follows from the simulations, is presented in Fig. 6. 
As follows from the figure, the flow reorganizes by diminishing PV gradients and 
flattening the radial PV distribution in the region between the circumpolar jet and 
the NPV , and thus tending to “cure” the instability, but this process is not yet 
over at t = 450/q” 1 . Indeed, opposite PV gradients of comparable strength are still 
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Figure 5. Snapshots of the pressure (left panel) and potential vorticity anomaly (right 
panel) at t = 150/o -1 during the saturation of the instability in the “jet-only” configu- 
ration. 



Figure 6. Left panel: Superposition of the observed mean azimuthal velocity profile with 
error margins (thin dashed), model azimuthal velocity profile (solid), and mean az- 
imuthal velocity profile (dashed) at t = 150/o -1 during the saturation of the instability 
in the “jet-only” configuration. Right panel: Evolution of the radial distribution of zon- 
ally averaged PVA in the “jet-only” configuration. 


present in the flow at this stage, thus maintaining the barotropic instability. We do 
not trace further the process of saturation which is well-known and leads to com- 
plete “curing” of the instability, appearance of counter-jets (although here we have 
a new element with respect to classical studies element, the 7- effect), and does 
not produce a quasi-stationary hexagonal structure. Note that the (diminishing) 
maximum of PVA shifts poleward during the saturation of the instability. 


Although the l = 6 instability is dominant, the growth rates of l = 5 and l = 7 
modes are rather close to the growth rate of l = 6 (Fig. 3, left bottom panel). If 
the initialization of direct numerical simulations is made with a combination of 
these modes with equal amplitudes and random phases, the velocity distribution 
is still meandering and possesses zonal counter-jets at the intermediate stages, and 
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Figure 7. Same as in the left panel of Fig 5, but for 
t = 350/ f 7 1 , (left panel), and 450 /q 1 (right panel). 

does not exhibit the hexagonal shape at late stages (not shown). 


4 Instability of the “jet + vortex” configuration and its nonlinear sat- 
uration 


It should be stressed that the “jet-only” approximation of section 3, which is 
being overwhelmingly used in the literature thus far, would be valid if the jet 
and the vortex were spatially well separated. As follows from Fig. 1, the velocity 
peaks are separated by approximately 3 Rd- A well-known property of the RSW 
is that the deformation radius Rd plays a role of a “screening” radius beyond 
which the velocity field generated by a vortex falls off exponentially. In spite of 
their rather well-separated velocity peaks, the velocity distributions of the central 
vortex and of the circumpolar jet are not sufficiently far from each other to prevent 
the central vortex influencing the perturbations of the jet. We will show below that 
this influence changes the evolution of the instability, although the most unstable 
mode of the composite jet + vortex system remains practically the same as in the 
“jet-only” case. 


4-1 Results of the linear stability analysis 


The most important result following from the linear stability analysis in the 
“jet+ vortex” configuration (Fig. 8) is that the most unstable mode has l = 6 
as in the “jet-only” case, and that its structure is practically the same. The big 
picture of the dependence of the azimuthal wavenumber of the most unstable mode 
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Re(co)= 0.044842, lm(co)= 0.01 7369 


Re(w)= 0.049225, lm((o)= 0.01 6855 




Figure 8. Structure of the most unstable mode with l = 6 of “jet-only ” configuration (left 
panel) and of “jet-fvortex” one (right panel)in the northern hemisphere corresponding 
to the background velocity profile of Fig.l. Vertical lines represent positions of critical 
levels. 



Figure 9. Same as in Fig. 4, but for the “jet-hvortex” configuration. The white square 
corresponds to the parameters in the right panel of Fig. 1 . 

on the parameter r 0 is also close to what was observed in the “jet-only” configu- 
ration, as follows from Fig. 10. The configuration corresponding to the South Pole 
has the most unstable mode with l = 8 (at the value Rd ~ 3200 km used as a 
reference for comparison with the “jet-only” case). The position of the southern 
polar jet being much farther from the pole in comparison with the northern polar 
jet, we can expect, as was the case in the “jet-only” configuration, that instability 
happens at higher l than the l = 6 leading to the hexagonal structure. The most 
unstable l can be pushed even higher by diminishing as follows from Fig. 11. 
The configuration corresponding to the South Pole has the most unstable mode 
with l = 8 (at the value Rd ~ 3200 km used as a reference for comparison with the 
“jet-only” case). The position of the southern polar jet being much farther from 
the pole in comparison with the northern polar jet, we can expect, as was the case 
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Figure 10. Same as in Fig. 4, but for the u jet+vortex v configuration. The white square 
corresponds to the parameters in the right panel of Fig. 1 . 





Figure 11. Linear growth rates o of different azimuthal modes of “jet+vortex” configu- 
ration at the North Pole (left panel) and South Pole (right panel) at different values for 
Rd. The non-dimensional velocity profiles are the same as Fig.l. 


in the “jet-only” configuration, that instability happens at higher l than the l = 6 
leading to the hexagonal structure. The most unstable l can be pushed even higher 
by diminishing as follows from Fig. 11. As one can infer from the observations 
(Antunano et ah, 2015) the southern circumpolar jet has a much higher azimuthal 
wavenumber than l = 8 (l = 10 — 11), we thus may conclude that the equivalent 
depth near the South Pole is, by some reason, smaller than in the vicinity of the 
North Pole. 
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4-2 Non-linear saturation of the instability 


Direct numerical simulations of the evolution of the instability of the “jet+vortex” 
system were performed using the same procedure as in the “jet-only” case. We 
present in Fig. 12 the snapshots of pressure and potential vorticity anomaly at 
times t = 150/o -1 and t = lOOO/o -1 clearly showing the formation and mainte- 
nance of a hexagonal structure matching the observed feature. Compared with the 
evolution of the “jet-only” configuration in Figs. 5, 7 the “jet+vortex” case ex- 
hibits much less deformation of the jet with time, resulting in the long-living quasi- 
stationary hexagonal pattern exhibiting sharp straight PV boundaries, which was 
not the case in the “jet-only” configuration. The evolution of the zonally averaged 
PVA is presented in Fig. 13. As follows from the figure, the jet+vortex configura- 
tion “cures” the instability at t — 1000 f^ 1 . Positive PV gradients become small 
and are shifted farther from the zone of negative PV gradients of the jet, which 
prevents from effective coupling of Rossby waves moving in opposite direction and, 
thus, eliminates the cause of instability. We performed a linear stability analysis of 
the mean zonal velocity profile corresponding to the PVA profile at t = 1000 ff 1 of 
Fig. 13 and found that although the unstable modes still exist, their growth rates 
are at least one order of magnitude less than those of unstable modes of initial 
profile. Hence, in accordance with qualitative considerations above, the instability 
does cure itself. This process, is, however, not sufficient to explain by itself why 
the hexagonal structure of finite amplitude persists. 

The inclusion of the NPV thus completely changes the evolution of the circumpolar 
jet. This could be further demonstrated by comparing the evolution of the mean 
azimuthal velocity in the “jet+vortex” case presented in Fig. 14 with that of the 
“jet-only” case in Fig. 6. The results are presented for two different intensities 
of the central vortex: one corresponding to the observations, and another one of 
the same shape but of higher amplitude. In the second case the system reaches a 
quasi-equilibrium configuration close to the observed one in high latitudes. Yet, a 
discrepancy is observed in low latitudes, where the velocity field in simulations does 
not match observations. It should be stressed that our goal is, in the first place, to 
understand the dynamical role of the polar vortex, so we did not seek to represent 
the details of the velocity field beyond the jet with our simple fit. There are also 
technical reasons for this decision. First, the 7 plane approximation is losing its 
accuracy while moving farther from the pole, so trying to reproduce fine details 
of the velocity profile does not make much sense. Second, in order to evacuate 
inertia- gravity waves, sponge boundary conditions were imposed at the boundary 
of the computational domain at r ~ 10. The mean velocity field, therefore should 
fall off to zero sufficiently far from the sponges, to be not perturbed by these latter. 
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Figure 12. Upper row: Snapshots of pressure distribution in “jet+vortex” configuration 
at time = 150, 1000/o -1 respectively, left and right panels. Lower row: corresponding 
potential vorticity anomaly. V(r) is the same as Fig.l. For better visibility the colorbars 
are not the same in two lower figures. 

However, the (weak) retrograde jet which is present in the data could play a role. 
For example, an exchange of angular momentum between retrograde and prograde 
jets in Saturn’s atmopshere was reported in (Liu and Schneider, 2015). Simulations 
on the whole sphere are needed in order to correctly reproduce this part of the 
velocity profile, which is beyond the scope of this paper. 

The same tendency is observed in the distribution of zonally averaged angular 
momentum represented in Fig. 15 by its deviation from the initial distribution. 
The angular momentum of the system is practically conserved (note that we have 
numerical dissipation and sponges at the boundaries). Fig. 15 gives an idea of 
angular momentum fluxes which are represented with the help of the difference of 
M{t) — M(to) the angular momentum density per unit mass M 

2 

M{r) = rv + f r — (4.1) 
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Figure 13. Evolution of the radial distribution of zonally averaged PVA in the 
“jet+vortex ” configuration. 




Figure 14. Superposition of the observed zonal velocity profile with error margins (thin 
dashed), the model zonal velocity profile used in initialization of numerical simulations 
(solid), and zonal velocity profile after the hexagon formation (dashed). In the right 
panel the amplitude of the initial vortex was taken deliberately higher (e = 0.2) than in 
the right panel (e = 0.15), in order to arrive to a close to observations end state. 

at a given time t and a fixed time to- We have chosen t 0 = 100 Jq 1 to avoid 
oscillations due to emission of inertia-gravity waves, and their partial reflection by 
(imperfect) sponges, at the initial stages of the evolution. The Figure shows a net 
angular momentum flux from the polar vortex to the jet, which diminishes in time 
in the jet+vortex configuration, but continues to act in the jet-only one. 


A tentative interpretation of the stabilization of the hexagonal structure of the 
circumpolar jet is that the central vortex is close enough to sweep the inner sec- 
ondary vortices appearing due to development of the barotropic instability of the 
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Time = 1 00 to 200 [f” 1 ] Time = 1 00 to 450 [f” 1 ] 



Figure 15. Difference of zonally averaged angular momentum density per unit mass M of 
“jet” and “jet-hvortex” configurations at times t = 200/q -1 (left panel), and t — 450/q" 1 
(right panel), and time t — 100 f^ 1 , in the northern hemisphere. Vertical dashed line 
represents the initial maximum of of the jet velocity. 



Figure 16. Evolution of the logarithms of the normalized amplitudes of the Fourier modes 
of the azimuthal velocity in a simulation of “jet+vortex” configuration in the northern 
hemisphere initialized with modes l = 2,3, ...,8 with the same (small) amplitudes and 
random phases. V(r) is the same as Fig.l. 


jet, and thus damp the meandering observed in the “jet-only” case. 

The above-described results correspond to initialization with a single dominant 
mode 1 = 6. We carried out simulations initialized with a spectrum of random- 
phase modes with different l = 2, 3,..., 8, which also leads to a hexagonal jet with 
similar properties as in the reference case of single-mode initialization, although 
it takes a longer time to settle. Fig. 16 illustrates the emergence of l = 6 as the 
dominant mode in this case. 
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5 Summary and discussion 


In the present paper we used the RSW model and the state-of-the-art Cassini ob- 
servations of the polar hexagonal jet and central vortex (Antunano et al., 2015) for 
a search of dynamical mechanisms that could explain the existence and longevity 
of the Saturn’s northern polar hexagon and the absence of its counterpart at the 
South pole. Insufficient knowledge of the vertical structure of the Saturn’s atmo- 
sphere justifies the use of such vertically integrated model, where all information 
about the vertical structure is condensed in a single parameter, the Rossby defor- 
mation radius R 4 . We used a fit of the observed velocity distributions to perform 
detailed linear stability analyses and to initialize nonlinear RSW simulations. In 
contrast with existing literature we analyzed not only a “jet-only”, but also a 
“jet+vortex” configuration, and showed the importance of including the central 
vortex in the system. 

Our main results are as follows 

(1) There indeed exists a region of parameters where the most unstable mode of 
the jet-only configuration has azimuthal wavenumber l = 6, although small 
changes of parameters may shift this value up or down. For the distance 
actually observed between the northern polar jet and the North pole, the 
azimuthal wavenumber l = 6 is obtained for a Rossby radius of deforma- 
tion Ra ~ 3200 km. 

(2) In the Southern hemisphere configuration, the most unstable mode has an 
azimuthal wavenumber higher than 8, with a precise value depending on Ra- 
This is related to the distance of the southern polar jet from the South pole 
being larger than in the Northern hemisphere. 

(3) The preceding results hold for both “jet-only” and “jet+vortex” cases. The 
difference between the two configurations is revealed by high-resolution fully 
nonlinear simulations initialized either with the most unstable mode of small 
amplitude, or with a random combination of unstable modes, superimposed 
onto the background velocity profile. In the “jet-only” configuration, nonlinear 
saturation of the instability leads to the formation of counter-jets and strong 
distortion of the hexagon at the late stages of the evolution, consistently 
with the standard scenario of saturation of the barotropic instability. On the 
contrary, in the “jet+vortex” configuration, the central vortex stabilizes the 
hexagon, leading at the late stages of the evolution to a configuration close 
to the observed one. 

The persistence of the hexagonal structure in our RSW simulations suggests that 
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Table 1 


Phase velocity of the hexagon 


e 

R d [km] 

c [r ad. /o], sole jet 

c [rad./o], jet+vortex 

0.08 

3200 

7.33 x 10 -3 

7.83 x 10 -3 

0.09 

3200 

8.40 x 10" 3 

8.88 x 10" 3 

0.10 

3200 

9.38 x 10 -3 

9.88 x 10“ 3 

0.11 

2400 

7.33 x 10 -3 

7.51 x 10“ 3 

0.12 

2400 

7.71 x 10" 3 

8.20 x 10" 3 

0.14 

2400 

9.06 x 10“ 3 

9.56 x 10“ 3 


the combination of central vortex and hexagonal jet structure could be an exact 
solution to the system. Proving this hypothesis is a challenge. Such proof would 
close the problem because, however robust and long-lived is this structure in our 
direct numerical simulations, although these latter are very long by the standards 
of computational fluid dynamics (more than a thousand of inertial periods), they 
remain far short of the time of observations of the Saturn’s hexagon. 

Compared to the existing literature, although the model and the background flow 
configurations are different, our “jet-only” configuration yields a “vortex-street” 
unstable northern polar hexagon similar to the results of Morales-Juberias et al. 
(2011), and the “jet+ vortex” configuration yields a “meandering jet” stable north- 
ern polar hexagon, similar to Morales-Juberias et al. (2015), featuring the sharp 
straight-line PV gradients on the six sides of the hexagon. This means that the pres- 
ence of the central vortex could stabilize the northern polar hexagonal jet already 
in a simple barotropic model, while a (baroclinic) vertical shear was necessary in 
Morales-Juberias et al. (2015) for stabilization of the hexagonal “meandering jet” . 


In spite of successful representation of the hexagon and explanation of its ab- 
sence at the opposite pole, we still have two caveats in our scenario. The first one 
is too high phase velocity of the hexagon. We present in Table 1 the values of 
the angular (or phase) velocity c of the hexagon for different values of parame- 
ters, as follows from the linear stability analysis of both the “jet-only” and the 
“jet+ vortex” configurations (it should be stressed that the angular velocities re- 
trieved from the direct nonlinear RSW simulations are close to results of the linear 
stability analysis). In both cases, the hexagon in our model is not quasi-stationary 
as in the observations (Sanchez-Lavega et al., 2014) (u = +0.0128 ± 0.0013 de- 
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gree longitude per Saturn day), but is rotating with c ~ +8 x 10 -3 radian. / 0 , 
i.e. uj = +6 degree longitude per Saturn day, meaning that the hexagonal pattern 
achieves a complete rotation in Saturn’s rotating frame in 60 days, which is much 
higher than in observations. Sanchez-Lavega et al. (2014) use, however, Voyager’s 
System III reference frame based on Saturn’s Kilometric Radiation (SKR). Not 
only Cassini’s SKR measurements yield a 8— min longer rotation period with re- 
spect to Voyager’s (Gurnett et ah, 2007), but alternative methods based on gravity 
measurements (Helled et al., 2015), or potential vorticity mapping (Read et ah, 
2009b), yield a 5- min to 7-min smaller rotation period, with respect to Voyager’s. 
Assuming our shallow-water simulations are consistent with the atmospheric dy- 
namics framework adopted in Read et al. (2009b), which entails that our simulated 
+6° /day drift rate would translate to a slower +3° /day drift rate in the Voyager 
SKR reference, and even a drift rate close to 0°/day in the Cassini SKR reference. 


The second one is the exaggerated zonal velocity in low latitudes resulting from 
the saturation of the instability. 

The second caveat is most probably due to the limitations of the gamma-plane 
model, simulations with the RSW model on the full sphere would be sufficient to 
answer the question whether this is a defect of the model, or of the gamma-plane 
approximation. The first one can be probably resolved by introducing a vertical 
structure in the model, such as an active lower layer, although this would require a 
number of ad-hoc hypotheses which we are reluctant to introduce. We should recall 
that baroclinic waves are slower than the barotropic ones, so a baroclinic structure 
would necessary slow down the system, the quantitative estimates depending on 
the details of the vertical structure. 

It should be also emphasized that we were working with a dissipationless version 
of the model. One of the reasons was to reduce the number of ad hoc hypotheses. 
Of course, an explicit (turbulent) viscosity and diffusion terms could be added in 
the r.h.s. of momentum and mass conservation equations of the model. In order to 
maintain any non-trivial steady state a forcing should be also added in this case 
to prevent dissipative decay. Both forcing and dissipation (the values of turbu- 
lent viscosity and diffusivity) would be necessarily ad hoc and should be adjusted 
to produce the observed mean flow. Above, we were pursuing the idea that the 
Saturn’s polar hexagon results from an instability of this mean flow, i. e. the po- 
lar jet plus polar vortex system. The only difference, if this system is considered 
as a solution of forced-dissipative equations instead of non-forced non-dissipative 
ones, is that eigen-frequencies of its small perturbations would be subject to ad- 
ditional viscous damping, which is necessarily weak in view of longevity of the 
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structure in question. A radiative damping can be also introduced in the system. 
In shallow- water modeling of the atmosphere such damping is usually represented 
by relaxation of the thickness of the layer towards some radiative equilibrium 
value with some characteristic relaxation time. As this time is presumably very 
long for Saturn, it was neglected. Another reason of working with inviscid model 
is that numerical scheme used for simulations of the saturation of the instability 
is practically dissipationless, which allows for very long-time runs. 

Radiative damping mentioned above is an example of diabatic flux in the system. 
Another one, which was recently invoked for explanation of the emergence of Sat- 
urn’s polar vortex with the help of the RSW model (O’Neill et al., 2016), and is 
recurrently evoked to explain the maintenance of alternating jets on giant planets, 
cf. e.g. (Showman, 2007) is moist convection. In relation to the discussion above, 
the moist waves are slower than the “dry” ones, so the moist effects could, in 
principle, explain the observed slowness of the hexagon’s rotation. The condensa- 
tion and moist convection can be included in the RSW model in a self-consistent 
way (Bouchut et al., 2009). For this, a new quantity, the bulk moisture content in 
the atmospheric column Q is introduced. It is conserved in the absence of phase 
transitions, but acquires a sink if condensation is switched on. The condensation 
happens whenever the threshold of saturation Q s is crossed. An extra convective 
flux through the upper surface of the atmospheric layer is added and linked to the 
latent heat of condensation. The resulting equations of so-called moist-convective 
rotating shallow water (mcRSW) read: 

/ 

d t v + (u*V)v + fz x v = —g'Vh, 

< d t h + V- (vh) = -f3P, (5.1) 

K dtQ + V-(Qv) = -P, 

with a coefficient /3 related to background stratification, and a simple relaxation 
parameterization for condensation P with relaxation time r: 

P= 9jl91 H(Q-Q s ). (5.2) 

T 

These new elements can be easily incorporated in the finite-volume numerical 
scheme we are using. 

Before making simulations with the full-strength mcRSW, we can get useful in- 
sights from analysing the evolution of Q in the absence of the condensation sink. 
Q then obeys the conservation equation, the third equation in (5.1) with P = 0. 
The results of previous RSW simulations with added passive tracer Q evolving 
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Figure 17. Snapshots of humidity initialized with uniform through the whole do- 
main value qo = 0.89 corresponding to the “jet+vortex” (left and middle panel at 
t — 150, 1000/o -1 ) and “jet-only” (right panel at t = 150/o _1 J configurations. 



Figure 18. Snapshots of condensation (thick black contours) and isobars 0.6, 0.7 and 0.75 
(grey contours ) in a most-precipitating environment, as follows from the simulations with 
moist-convective RSW of Bouchut et a 1. (2009). Humidity is distributed uniformly in the 
Bow domain, with initial value Q = 0.48 with a saturation threshold at Q s — 0.51. Other 
initial or saturation values do not change the evolution qualitatively. 

from initially uniform distribution are presented in Fig. 17, for both “jet-only” 
and “jet+vortex” configurations. Fig. 17 shows that the excess of humidity re- 
sulting from the evolution of the system is concentrated in the vicinity of the jet, 
and does not affect the central vortex. It is generally known, and confirmed by 
the simulations with moist-convective RSW in Lambaerts et al. (2011); Rostami 
and Zeitlin (2017), that condensation and related latent heat release lead to in- 
tensification of the cyclones and damping of the anticyclones So we can infer that 
in the absence of additional sources of moisture, condensation would not help to 
maintain the central vortex, and thus to increase the longevity of the system, but 
would only impact the hexagonal jet on its equatorward flank. Simulations with 
the full moist-convective RSW (5.1) confirm this prediction, as follows from Fig. 
18 where evolution of condensation in time is displayed. Again, this result may 
change in the two-layer version of the model, but its introduction requires more 
information on the vertical structure of the Saturn’s atmosphere. 
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